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Abstract 

Short description of the FORTRAN-codes for an analysis of the ul- 
trashort pulse dynamics is presented. We consider: 1) the aberration-less 
approximation and the momentum method for the search of the single 
pulse stability regions in the laser with the soft-aperture Kerr-lens mode 
locking; 2) the distributed complex Ginzburg-Landau model for the same 
aim; 3) the generalized Schrodinger equation for the analysis of the fem- 
tosecond pulse propagation in the tapered and photonic crystal fibers. 

1 Introduction 

Analysis of the intra- and extra-cavity pulse dynamics has the common features: 
it can be based on the distributed models of the complex Ginzburg-Landau type 
P or the generalized Schrodinger type [2]. The basic difference of ultrashort 
pulse laser dynamics from, for example, intra-fiber propagation of the femtosec- 
ond pulses is the essential contribution of the dissipation into former. However, 
for the direct simulation the appearance of the non-Hamiltonian part in the 
dynamical equation doesn't cause some problems. It allows considering the dis- 
tributed models on the common basis (some physical and formal backgrounds 
can be found in 0). 

The main problem results from the transition to the higher-dimensional mod- 
els. The trivial dimension is 1 + 1, i.e. propagation distance (or cavity transit 
number for a laser) + local time (for an analysis of the ultrashort pulse dy- 
namics). Our experience proved that such dimension is quite sufficient for the 
solution of the various problems, for example, spectral continuum generation 
in the fibers or multiple pulse operation of the femtosecond lasers. This ap- 
proach adequately describes the pulse dynamics and at the same time is not 
too complicated to remain physically meaningful. The minimum set of the gov- 
erning parameters allows the overall optimization by the direct scanning of the 
parametrical space 0] E] • 

However, a consideration of the self-start ability of the soft-aperture Kerr- 
lens mode- locked lasers (see [2]) and an analysis of the real- world laser config- 



urations requires at the least 1 + 2 dimensions, where one transverse spatial 
dimension appears as a result of the rotational symmetry. The direct simu- 
lations on this way exceed the resources of the desk-top computer, which is 
common tool for the laser community. Therefore we have to reduce the dimen- 
sion, for example, by the means of the so-called aberration-less approximation. 
This allows performing the simulation not on the time-spatial grid but on the 
grid formed by the parametrical set defining the trial (aberration-less) solution. 

Below we consider some examples of 1 + 1 and 1 + 2 dimensional models. 
The latter model has in fact reduced dimension due to the use of the momentum 
method. It should be noted that the analytical computations underlying the 
numerical codes were realized in the MAPLE computer algebra system 

2 Passive mode-locking and nonlinear complex 
Ginzburg-Landau equation 

This approach is based on 1+1 dimensional model in the framework of the 
so-called nonlinear Ginzburg-Landau equation, which describes the Kerr-lens 
mode locking as an action of the fast saturable absorber governed by the few 
physically meaningful parameters, viz., its modulation depth 7 and the inverse 
saturation intensity a. 

The master equation describing the ultrashort pulse generation in the Kerr- 
lens mode-locked solid-state laser is: 
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where a(z, t) is the field amplitude (so that |a| 2 has a dimension of the intensity), 
z is the longitudinal coordinate normalized to the cavity length, t is the local 
time, a is the saturated gain coefficient, p is the linear net-loss coefficient tak- 
ing into account the intracavity and output losses, tf is the group delay caused 
by the spectral filtering within the cavity, (3 m are the m-order group-delay dis- 
persion coefficients, S = Igniw^jc = 2Tm2l g /(\on) is the self-phase modulation 
coefficient, uiq and Ao are the frequency and wavelength corresponding to the 
minimum spectral loss, n and rii are the linear and nonlinear refraction coeffi- 
cients, respectively, l g is the double length of the gain medium (we suppose that 
the gain medium gives a main contribution to the self-phase modulation). The 
last term in Eq. (1) describes the self-steepening effect and for the simplification 
will be not taken into account in the simulations. As an additional simplification 
we neglect the stimulated Raman scattering in the active medium 0. 

For the numerical simulations in the framework of the distributed model it 
is convenient to normalize the time and the intensity to tf = Aq/(AAc) and 1/(5, 
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respectively ( AA is the gain bandwidth) . The simulation were performed on the 
2 12 x 10 4 mesh. Only steady-state pulses were considered. As the criterion of 
the steady-state operation we chose the peak intensity change less than 1% over 
last 1000 cavity transits. 

Note that the local time interval, which is equal to the cavity period w 10 
nanoseconds, is not covered in our case (2 12 x t f w 20 -7- 100 picoseconds). 
This puts the questions about stability against the multipulsing with the large 
inter-pulse separations. Additionally we can not be sure in the ability of the 
spontaneous appearance of the mode locking (problem of the mode locking self- 
start ability). 

The solution of Eq. (1) is based on the fast Fourier-transform split-step 
method (see Appendix 1). We symmetrized non-Hamiltonian (square brackets 
in Eq. (1)) and Hamiltonian part (braces in Eq. (1)) separately. The mode 
locking in the considered model is governed by the only four basic parameters: 
a — Pi 02, 7, and a. This allows unambiguous multiparametric optimization. 
In the presence of the higher-order dispersions, the additional (3 m parameters 
appear. This complicates the optimization procedure, but keeps its physical 
clarity. As an initial condition we take the analytical solution of the cubic 
Ginzburg-Landau or Schrodinger equation 3 . 

Some results obtained on the basis of this model are presented in 0] |SJ [7] . 

3 Spectral continuum generation in the tapered 
fiber 

Generation of spectral supercontinuum became a hot topic in optics in recent 
years [8]. In the crystal or tapered fiber the propagating field has a compar- 
atively determinate spatial structure due to a strong confinement. Therefore, 
the 1 + 1-dimensional simulations give a quite thorough result. However, the 
nonlinearity in such fibers is enhanced by their small core size. As a result, a 
set of the nonlinearities is wider than that described by the Hamiltonian part 
of Eq. (1), which is the high-order nonlinear Schrodinger equation. The non- 
trivial generalization can be obtain due to taking into account the stimulated 
Raman scattering j^j: 

dz ^ ml dt m v ' 

m>2 

/ 

-Sfa J R{t)\a(z,t-t)\ 2 dt, 

— oo 

where S — uhoq/c is the self-phase modulation coefficient, (3 m is the mth- 
order group- velocity dispersion coefficient, / is the fraction of the stimulated 
Raman scattering contribution to the nonlinear refractive index of the fiber, 
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R(t) = ^r^j exp (^—^f-j sin yj^j is the Raman response function |1U1 IX 1| . 
T\ = 12.2 fs and T2 = 32 fs define the phonon oscillation period and its dumping 
time, respectively. 

We normalized t to 1 fs (the normalization to the initial pulse width is 
convenient, too), z to the nonlinear length L n i — (ujq^Io/c) defined by the 
initial pulse intensity Iq . The simulations were carried out on the mesh with the 
time step 1 fs (2 13 points) and the spatial step 10 -3 L n \. The solution of Eq. (2) 
was based on the fast Fourier-transform split-step method with the evaluation 
of the Raman response in the time domain (see Appendix 2) . 

The analysis of the pulse propagation in the tapered fibers requires the 
attention to the so-called transient sectors: before (after) the tapered sector 
with the almost constant waist there is the convergent (divergent) sector, where 
the fiber characteristics change from those in the single-mode fiber to those in 
the tapered fiber (or vice versa). The exact law of these changes is unknown, 
but we used the simple linear approximation for the evolution of the intensity 
and the dispersion coefficients. Note, that the normalization of the intensity 
was defined through the parameters of the tapered sector. 

4 Aberration-less approximation: analysis of the 
real-world laser configurations 

Above the problem of the ultrashort pulse stability in the mode-locked laser 
were formulated on the basis of the distributed 1 + 1-dimensional Ginzburg- 
Landau model. We noted also some problems of such model. Here we present 
an analysis on the basis of the time-spatial model. The spatial distribution for 
the laser beam is assumed to be Gaussian that reduces the problem to 1 + 2- 
dimensions. The free-space propagation of the Gaussian beam can be considered 
on the basis of the usual ABCD-matrix formalism while the propagation 
inside the nonlinear active medium is described by the following equation: 
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Here (3' 2 and t'^ are the group-velocity dispersion and the inverse group-velocity 
delay coefficients (for ZnSe laser we used /3 2 2 =2054 fs 2 /cm and £^=13 fs/cm 
0]). The left-hand side of Eq. (3) describes the non-dissipative factors: thermo- 
lcnsing (i9—k^-(P a exp (—(z)/ (innaKth), k is the wave number, ^f- is the co- 
efficient of the refractive index thermo-dependence (5.35x lO^K^ 1 for ZnSe), 
£ is the loss coefficient at the pump wavelength, P a is the pump power, n t h is 
the thermo-conductivity coefficient (0.172 WK^cmT 1 for ZnSe)); diffraction 
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(in the cylindrically symmetrical case) ; group- velocity dispersion and self-phase 
modulation (providing self-focusing for radially varying beam, x^n^kjriQ). The 
right-hand side of Eq. (3) describes the dissipative factors inside the gain 
medium: radially varying gain (providing gain guiding and soft aperture ac- 
tion, a and w p are the saturated gain coefficient and the pump beam size, 
respectively); spectral filtering caused by the gain band profile. The saturated 
gain can be expressed in the following way: 

a = 7 r, (4) 

fa nw 2 ( + + A 

where v=Eirw 2 / (2P g T p ) (P g is the generation power, E is the generation energy, 
uj p is the pump frequency, r p is the pulse width, T cav is the cavity period, a a is 
the absorption cross-section, I s is the gain saturation energy), w is the genera- 
tion mode beam size, v = \Jix /I for the pulse with the Gaussian time-profile, 2 
for the sec/i-shaped pulse and 1 for the CW (in the latter case T p =T cav ). The ap- 
proximated solution of Eq. (3) is based on the so-called aberration-less approx- 
imation: the propagating field has the invariable spatial-time profile, which is 
described by the set of the z-dependent parameters. In the non-dissipative case 
this approximation allows the variational approach providing rigorous descrip- 
tion of the Gaussian beam propagation outside the parabolical approximation 

In the dissipative case we use the momentum method 14 and consider 
the momentums resulting from the symmetries of Eq. (3). The a — > aexp {i<j>) 
invariance, the transverse and time translating invariance suggest the following 
momentums |15j : 



T m ,n= I I r m t n \a\ 2 drdt, (5) 

,„ ,„ , da* . On \ 

r t I a— a — drdt. 

or or ) 

Like the variational approach we can substitute to Eqs. (3, 5) the trial so- 
lution describing the ultrashort pulse. If we take the Gaussian time-spatial 

profile a (z, r, t) = W(r) exp ^(r)-^^ + ib(r)r 2 - ^ + i^{r)t 2 ) (W(r) 

is the complex amplitude, 2w' 2 =w 2 , G(r) is the pulse amplification parameter 
excepting the geometrical magnification for the Gaussian beam) , the equations 
describing the evolution of the pulse and beam parameters are |16| : 
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where P$ and w' are the power and the beam size before the active medium, 
respectively. This system can be solved on the basis of the fourth-order Runge- 
Kutta method (see Appendix 3; beta, alpha, gamma, delta and psi correspond 
to the right-hand side of Eqs. (6) for b, w, G, t and ip, respectively). 

5 Conclusion 

Above considered models allow the different generalizations. For example, in 
the framework of the Ginzburg-Landau model the stimulated Raman scattering 
inside the active medium can be taken into account immediately (by analogy 
with the fiber optics). The high-order nonlinear Schrodinger equation can be 
generalized in order to take into account the birefringence. The momentum 
method requires an additional analysis for the dissipative propagation. Nev- 
ertheless, in the presented form it can be used for the Kerr-lens mode-locked 
lasers optimization. Outside the aberration-less approximation for the field 
time-profile, the model allows most adequate description of the mode-locked 
lasers, but the computational time can be enormous in this case (for up-to-date 
desk-top computers). 

The considered numerical codes were prepared as a result of the prelimi- 
nary analysis in the framework of the computer algebra system MAPLE (see 
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Appendix 1: 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c ULTRASHORT PULSE STABILITY 

c IN THE KERR-LENS MODE-LOCKED LASER: 

c ANALYSIS ON THE BASIS OF THE COMPLEX 

c GINZBURG-LANDAU EQUATION 

c 

c V.L.Kalashnikov 

c Photonics Institute, Technical University of Vienna 
c e-mail: kalashnikov@tuwien.ac. at 
c web-site: 



http://www.geocities.com/optomaplev 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

c PARAMETERS: 
c 

c alpha is the saturated net-gain; 

c beta is the group-delay dispersion; 

c D3 is the third-order dispersion; 

c gamma is the fast absorber modulation depth; 

c delta is the fast absorber saturation parameter. 

c Vr and Vi are the real and imaginary parts of the field, respectively; 

c V is the field intensity; VV is the spectral intensity; 

c Vmax is the maximum pulse intensity; 

c En is the generation energy; 

c width is the pulse width; 

c shift is the spectrum maximum shift. 

c 

c The number of the considered steady-state pulses is defined by 

c the number counter 

c 

c All values are dimensionless: the intensity is normalized to the self- 

c phase modulation coefficient, the time is normalized to the inverse 

c bandwidth of the spectral filter, the propagation distance is normalized 

c to the cavity length 

REAL*8 Vr(4096),Vi(4096),V(4096),VV(4096),alpha,beta,gamma,delta 
common /comin/ts(4096),ntab 

DOUBLE PRECISION rho,rhol,rho2,tau,C,S,X,Y,Argum,Fout,0mega,Vm 
DATA Nt,Nst,Pi/4096,12,3.14159265358979323846264338d0/ 



epsilon = l.e-4 
open(l,file= 'Landau_G.dat') 
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write(l,*)'gamma, alpha, beta, D3, delta, Vmax, En, width, shift' 
gamma = 0.05 
D3 = -150. 

do j2=l,10 

alpha= (gamma-epsilon) * j 2 / 1 OdO 
do j3=l,100 
bcta=-j3 
do j4=l,50 

delta=10**(-2dO+4dO*j4/50dO) 

Fout = IdO/Nt 

c Initialization of the initial field 

K = 

tau=sqrt (gamma-alpha) 

rhol = sqrt(2d0)*tau/sqrt(gamma*delta) !from cubic Landau- Ginzburg 

rho2 = sqrt(-beta)*tau ! from Schrodinger 

if (r ho 1 . gt . r ho 2 ) t hen 

rho=rhol 

else 

rho=rho2 
end if 

do j=l,Nt 

Vr(j)=rho/cosh((j-2048)*tau) 

Vi(j)=0. 
end do 

1 K = K+l 

c DISSIPATIVE PART 

c First amplification step 
do I=l,Nt 

Vr(i) = Vr(i)*exp(0.5d0*alpha) 
Vi(i) = Vi(i)*exp(0.5d0*alpha) 
end do 
c First filter action 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

Argum=Omega**2 ! gain profile and action 
Vr(I)=Vr(I)*exp(-0.5dO*Argum) 
Vi(I)=Vi(I)*exp(-0.5dO*Argum) 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
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First nonlinear part's action 
do I=l,Nt 

Argum=gamma/(ldO + delta*(Vr(I)**2+Vi(I)**2)) 
Vr(I)=Vr(I)*exp(-0.5dO*Argum) 
Vi(I)=Vi(I)*exp(-0.5dO*Argum) 

end do 
Second amplification step 
do I=l,Nt 

Vr(i) = Vr(i)*cxp(0.5d0*alpha) 
Vi(i) = Vi(i)*exp(0.5d0*alpha) 
end do 
c Second filter action 

call fftmn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(LLE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omcga=2.*Pi*Is*Fout 

Argum=Omcga**2 ! gain profile and action 
Vr(I)=Vr(I)*exp(-0.5dO*Argum) 
Vi(I)=Vi(I)*exp(-0.5dO*Argum) 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c Second nonlinear part's action 
do I=l,Nt 

Argum=gamma/(ldO + dclta*(Vr(I)**2+Vi(I)**2)) 
Vr(I)=Vr(I)*exp(-0.5dO*Argum) 
Vi(I)=Vi(I)*exp(-0.5dO*Argum) 

end do 
c HAMILTONIAN PART 
c First dispersion step 

call fftmn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(0.5dO*((beta/2dO)*Omega* s,! 2 - (D3/6dO)*Omega**3)) 
S=sm(0.5dO*((beta/2dO)*Omega* s,! 2 - (D3/6dO)*Omega**3)) 
X=Vr(I) 
Y=Vi(I) 

Vr(I)=X*C+Y*S 
Vi(I)=Y*C-X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c First self-phase modulation step 
do I=l,Nt 

9 



c 



c 



Argum=0.5dO*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C+Y*S 
Vi(I)=Y*C-X*S 
end do 

c Second dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(0.5dO*((beta/2dO)*Omega**2 - (D3/6dO)*Omcga**3)) 
S=sin(0.5dO*((beta/2dO)*Omega**2 - (D3/6dO)*Omega**3)) 
X=Vr(I) 
Y=Vi(I) 

Vr(I)=X*C+Y*S 
Vi(I)=Y*C-X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c Second self-phase modulation step 
do I=l,Nt 

Argum=0.5dO*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C+Y*S 
Vi(I)=Y*C-X*S 
end do 

c Analysis of Generation Field 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 

do i=l,nt 

if(i.lc.nt/2)then 

j=nt/2+i 

VV(j)=Vr(i)**2 + Vi(i)**2 
else 

j=i-nt/2 

VV(j)=Vr(i)**2 + Vi(i)**2 
end if 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 

Sm=W(l) 

do 6 1=1, NT 
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if(W(I).LT.Sm)goto 6 
Sm=VV(I) 
Ism=I 
6 continue 

shift=Ism-Nt/2 
do I=1,NT 

V(I)=Vr(I)**2+Vi(I)**2 
end do 



Vm=V(l) 

do 2 I=1,NT 

if(V(I).LT.Vm)goto 2 

Vm=V(I) 

Im=I 
2 continue 

if(k.eq.9000)cont_p=Vm 

if(Vm.lt.le-10)goto 5 

if(k.lt.l0000)goto 1 
c Analysis of Generation Field 

do I=1,NT 

IF(V(I).ge.Vm/2.and.V(I-l).le.Vm/2)hl=(2.*I-l)/2. 
IF(V(I).le.Vm/2.and.V(I-l).gc.Vm/2)h2=(2.*I-l)/2. 
end do 

En = 0. 

do I=l,Nt 

En=En+V(i) 

end do 

if(Vm.gt.0.)then 
stab=abs(cont_p-Vm) /Vm 
else 

stab=l. 

end if 
Vmax=Vm 
c PULSE NUMBER COUNTER 
call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 

c Filter 

do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 
IF(I.GE.Nt/2+2)Is=I-l-Nt 
Omega=2.*Pi*Is*Fout 
Argum= (0. 1 * (h2-hl ) ) * * 2* Omega* *2 
Vr(I)=Vr(I)*exp(-Argum) 
Vi(I) = Vi(I) *cxp (- Argum) 
end do 
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call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 



do I=1,NT 

V(I)=Vr(I)**2+Vi(I)**2 

end do 

Vm=V(l) 

do 3 1=1, NT 

if(V(I).LT.Vm)goto 3 

Vm=V(I) 

Im=I 

3 continue 

nm=0 

do i=2,nt-l 

if(V(i).ge.V(i-l).and.V(i).gt.Vm/10.)then 

if(V(i+l).lt.V(i))nm=nm+l 

else 

end if 

end do 

width=(h2-hl) 
if(nm.eq.l.and.stab.lt.0.01)then 

writc(l,fmt=4)gamma, alpha, beta, D3, delta, Vmax, En, width, shift 

4 format(f4.2,lx,f5.3,lx,f6.1,lx,f6.1,lx,f5.1,lx,f7.4,lx,f7.3,lx, 
#f7.1,lx,f7.1) 

else 
end if 

5 end do 

end do 
end do 

close(l) 
end 

c Fast Fourier transformation 

SUBROUTINE FFTINN(XB,XM,NT,NST,SIGN) 
REAL*8 XB(NT),XM(NT) 
COMMON/COMIN/TS(4096),NTAB 

DOUBLE PRECISION XB1,XB2,XM1,XM2,PRIR,W1,W2,FP,TS 

IF(NTAB.GT.0)GOTO 30 

NTAB=4096 

PRIR=3.14159265358979323846264338d0/NTAB 
DO 20 I=1,NTAB 
FP=PRIR*(I-1) 
20 TS(I)=SIN(FP) 
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30 CONTINUE 
LEC1=NT-1 
LEC2=NT/2 
J=l 

DO 10 I=1,LEC1 

IF(I.GE.J)GO TO 8 

XB1=XB(J) 

XM1=XM(J) 

XB(J)=XB(I) 

XM(J)=XM(I) 

XB(I)=XB1 

XM(I)=XM1 

8 L=LEC2 

9 IF(L.GE.J)GOTO 10 
J=J-L 

L=L/2 
GO TO 9 

10 J=J+L 
INCP=NTAB*2 
NSDV=NTAB/2 
JLI=1 

JKI=1 

KLI=NT 

DO 3 I=1,NST 

JKI=JKI+JKI 

KLI=KLI/2 

INCP=INCP/2 

INI=1 

DO 2 J=1,JLI 
LEC1=J 
W2=TS(INI) 
IF(INI-NSDV)5,5,6 

5 W1=TS(INI+NSDV) 
GO TO 7 

6 W1=-TS(INI-NSDV) 

7 IF(SIGN.GT.0.)W2=-W2 
DO 1 K=1,KLI 
LEC2=LEC1+JLI 
XB1=XB(LEC1) 
XM1=XM(LEC1) 

XB2=W1*XB(LEC2)-W2*XM(LEC2) 

XM2=W1*XM(LEC2)+W2*XB(LEC2) 

XB(LEC1)=XB1+XB2 

XM(LEC1)=XM1+XM2 

XB(LEC2)=XB1-XB2 

XM(LEC2)=XM1-XM2 
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1 LEC1=LEC1+JKI 

2 INI=INI+INCP 

3 JLI=JLI+JLI 
IF(SIGN.LT.O.)RETURN 
FP=1./NT 

DO 4 I=1,NT 
XB(I)=XB(I)*FP 

4 XM(I)=XM(I)*FP 
RETURN 

END 



Appendix 2: 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cSPECTRAL CONTINUUM GENERATION IN THE TAPERED FIBER: 
c GENERALIZED SCHROEDINGER EQUATION 

c 

c V.L.Kalashnikov 
c Photonics Institute, Technical University of Vienna 

c e-mail: kalashnikov@tuwien.ac.at 

c web-site: http://www.geocities.com/optomaplev 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

c The description of the program can be found on http: / /www. geocities.com/ 
c optomaplev/programs/fortran.html 

c 

c PARAMETERS: 

c 

c Amp is the initial pulse intensity (1 for the tapered section) 
c ar_begin is defined by the ration of the tapered fiber and single-mode 
c fiber cross-sections 
c width is the initial pulse width (in fs) 

c D_begin and D3_begin are the group-velocity and third-order dispersions 
c of the single-mode fiber 

c D_end and D3_end are the group-velocity and third-order dispersions 
c of the tapered fiber 

c All dispersions are defined through the dispersion lengths normalized to 
c the nonlinear length in the tapered fiber 

c N_steps and N_t are 1000*L/Lnl, where L is the length of the transient 
c fiber sector or the tapered sector, respectively Lnl is the 
c nonlinear length of the single-mode or tipered fiber, respectively 

REAL*8 Vr(8192),Vi(8192),V(8192),R(8192) 
REAL*8 Vr.out(8192),VLout(8192) 
INTEGER N.steps,N_t 
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common /comin/ts(8192),ntab 



DOUBLE PRECISION Amp,width,C,S,X,Y,Argum,Fout,Omega,Vm,n2eff 
DATA Nt,Nst,Pi,Dx/8192,13,3.14159265358979323846264338d0,ld-3/ 
DATA Pr,Tl,T2/0.15d0,12.2d0,32d0/ 



open(l,file='phase.dat') 
opcn(2,filc= 'spcctrum.dat') 
open(3,file='intensity.dat') 
Fout = IdO/Nt 
c Raman responce function 
do i=l,Nt 

R(i) = ((Tl**2+T2**2)/(Tl*T2**2))*cxp(-(i-l)/T2)*sin((i-l)/Tl) 
end do 

c Initialization of the initial field 

K = 

width=28.36d0 



D.begin = -width**2/7.55d0 
D.end = width**2/0.53d0 
D3.begin = -0.58dO*D_begin*width 
D3.cnd = -0.78dO*D_end*width 
ar.begin = 0.0149 
N_steps = 3774 
N_t = 16992 

AAA = ar_begin**(-l./(2.*Nj3teps)) 

D = D.begin 
D3 = D3_begin 

Amp = sqrt(ar .begin) 

do j=l,Nt 

Vr(j)=Amp/cosh((j-4096)/width) 

Vi(j)=0. 
end do 
do j=l,Nt 

write(3,*)k,j,Vr(j)**2+Vi(j)**2 
end do 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do i=l,nt 
if(i.le.nt/2)thcn 
j=nt/2+i 

V(j)=Vr(i)**2 + Vi(i)**2 
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Vr_out(j)=Vr(i) 
VLout(j)=Vi(i) 
else 

j=i-nt/2 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

Vi_out(j)=Vi(i) 

end if 

end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
do i=l,Nt 

write(l,*)k,i,Vr_out(i),Vi_out(i) 
end do 
do j=l,Nt 
write(2,*)k,j,V(j) 
end do 

c FIRST TRANSIENT SECTOR 
1 K = K+l 

writc(*,*)K 

c First dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omcga=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega* 
S=sin(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega*' 
X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c First self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO*(ldO-Fr)*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c First Raman step 
do I=l,Nt 
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Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J)**2) 
end do 

C=cos(Dx*0.5dO*Fr* Argum) 

S=sin(Dx*0.5dO*Fr*Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 

end do 

c Second dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omcga=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omcga**2 - (D3/6dO)*Omega* 
S=sin(Dx*0.5dO*((D/2dO)*Omcga**2 - (D3/6d0) * Omega* >: 
X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c Second self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO*(ldO-Fr)*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c Second Raman step 
do I=l,Nt 

Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J)**2) 
end do 

C=cos(Dx*0.5dO*Fr*Argum) 

S=sin(Dx*0.5dO*Fr*Argum) 
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X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
do I=l,Nt 
Vr(i) = Vr(i)*AAA 
Vi(i) = Vi(i)*AAA 
end do 

D = D.bcgin + K*(D_cnd - D_bcgin)/N_steps 
D3 = D3_bcgin + K*(D3_end - D3_begin)/N_steps 
if ( (K / 1000) *1000.eq.K)thcn 
do j=l,Nt 

write(3,*)k,j,Vr(j)**2+Vi(j)**2 
end do 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do i=l,nt 
if(i.le.nt/2)then 
j=nt/2+i 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

Vi_out(j)=Vi(i) 

else 

j=i-nt/2 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

Vi.out(j)=Vi(i) 

end if 

end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
do i=l,Nt 

write(l,*)k,i,Vr_out(i),Vi_out(i) 

end do 

do j=l,Nt 

write(2,*)k,j,V(j) 

end do 

else 

end if 

if(k.lt.N_steps)goto 1 
c TAPERED SECTOR 

K = K+l 

writc(*,*)K 

D = D_cnd 
D3 = D3_end 

c First dispersion step 



18 



call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega**3)) 
S=sin(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega**3)) 
X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c First self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO*(ldO-Fr)*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S— sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c First Raman step 
do I=l,Nt 

Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J)**2) 
end do 

C=cos(Dx*0.5dO*Fr*Argum) 

S=sin(Dx*0.5dO*Fr*Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 

end do 

c Second dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omcga**2 - (D3/6dO)*Omega**3)) 
S=sin(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega**3)) 
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X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c Second self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO* (IdO-Fr) * (Vr (I)**2+Vi(I) * 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c Second Raman step 
do I=l,Nt 

Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J) 
end do 

C=cos(Dx*0.5dO*Fr*Argum) 
S=sin(Dx*0.5dO*Fr* Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 

if((K/1000)*1000.eq.K)thcn 
do j=l,Nt 

write(3,*)k,j,Vr(j)**2+Vi(j)**2 
end do 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do i=l,nt 
if(i.lc.nt/2)thcn 
j=nt/2+i 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

Vi_out(j)=Vi(i) 

else 

j=i-nt/2 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

VLout(j)=Vi(i) 
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end if 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
do i=l,Nt 

write(l,*)k,i,Vr_out(i),Vi_out(i) 

end do 

do j=l,Nt 

write(2,*)k,j,V(j) 

end do 

else 

end if 

if(k.lt.N_steps+N_t)goto 2 
c SECOND TRANSIENT SECTOR 

K = K+l 

writc(*,*)K 
c First dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega* 
S=sin(Dx*0.5dO*((D/2dO)*Omcga**2 - (D3/6dO)*Omega*^ 
X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c First self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO*(ldO-Fr)*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c First Raman step 
do I=l,Nt 

Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J)**2) 
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end do 

C=cos(Dx*0.5dO*Fr*Argum) 

S=sin(Dx*0.5dO*Fr*Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 

end do 

c Second dispersion step 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do I=l,Nt 

IF(I.LE.Nt/2+l)Is=I-l 

IF(I.GE.Nt/2+2)Is=I-l-Nt 

Omega=2.*Pi*Is*Fout 

C=cos(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omega**3)) 
S=sin(Dx*0.5dO*((D/2dO)*Omega**2 - (D3/6dO)*Omcga**3)) 
X=Vr(I) 
Y=Vi(I) 

Vr(I)= X*C - Y*S 
Vi(I)= Y*C + X*S 
end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
c Second self-phase modulation step 
do I=l,Nt 

Argum=Dx*0.5dO*(ldO-Fr)*(Vr(I)**2+Vi(I)**2) 
C=cos(Argum) 

S=sin(Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
end do 
c Second Raman step 
do I=l,Nt 

Argum = OdO 
do J=1,I 

Argum=Argum + R(I-J+l)*(Vr(J)**2+Vi(J)**2) 
end do 

C=cos(Dx*0.5dO*Fr*Argum) 

S=sin(Dx*0.5dO*Fr*Argum) 

X=Vr(I) 
Y=Vi(I) 
Vr(I)=X*C - Y*S 
Vi(I)=Y*C + X*S 
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end do 

do I=l,Nt 

Vr(i) = Vr(i)/AAA 

Vi(i) = Vi(i)/AAA 

end do 

D = D.cnd + (K-N.steps-N_t)*(D_begin - D.cnd)/N_steps 
D3 = D3_end + (K-N_steps-N_t)*(D3_begin - D3_end)/N_steps 
if((K/1000)*1000.eq.K)then 
do j=l,Nt 

write(3,*)kj,Vr(j)**2+Vi(j)**2 
end do 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do i=l,nt 
if(i.lc.nt/2)thcn 
j=nt/2+i 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

VLout(j)=Vi(i) 

else 

j=i-nt/2 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

Vi_out(j)=Vi(i) 

end if 

end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
do i=l,Nt 

write(l,*)k,i,Vr_out(i),VLout(i) 

end do 

do j=l,Nt 

write(2,*)k,j,V(j) 

end do 

else 

end if 

if(k.lt.N_t + 2*N_steps)goto 3 
do j=l,Nt 

writc(3,*)k,j,Vr(j)**2+Vi(j)**2 
end do 

call fftinn(Vr,Vi,Nt,Nst,l.) ! from Time to Frequency 
do i=l,nt 
if(i.le.nt/2)then 
j=nt/2+i 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

VLout(j)=Vi(i) 

else 
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j=i-nt/2 

V(j)=Vr(i)**2 + Vi(i)**2 

Vr_out(j)=Vr(i) 

VLout(j)=Vi(i) 

end if 

end do 

call fftinn(Vr,Vi,Nt,Nst,-l.) ! from Frequency to Time 
do i=l,Nt 

write(l,*)k,i,Vr_out(i),Vi_out(i) 
end do 
do j=l,Nt 
write(2*)k,j,V(j) 
end do 

close(l) 
close(2) 
close(3) 
end 

c Fast Fourier transformation 
see Appendix 1 



9 Appendix 3 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c ULTRASHORT PULSE STABILITY IN THE KERR-LENS 
c MODE-LOCKED LASER: ANALYSIS ON THE BASIS OF 
c THE MOMENTUM METHOD 
c 

c V.L.Kalashnikov 

c Photonics Institute, Technical University of Vienna 

c e-mail: kalashnikov@tuwien.ac.at 

c web-site: http://www.geocities.com/optomaplev 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c This program is realized on the basis of the computer algebra approach 
c (the corresponding description can be found on http:/ /www. geocities.com/ 
c optomaplev/programs/fortran.html) 
c 

c PARAMETERS: 
c 

c a is the distance from the out-put plane mirror to the first folding 
c mirror; 

c c is the distance from the totally reflecting plane mirror to the 

c second folding mirror; 

c b is the distance between folding mirrors; 
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bs and bf give the limits of the scanning on b; 

bl is the distance of the active medium facet from the first folding 

mirror; 

bis is the starting bl 

f is the focus length of the folding mirrors; 
z is the active medium length; 
n is the refractive index of the active medium; 
lambda is the generation wavelength; 

epsl gives the criterion of the convergence to the steady-state 
solution; 

P is the pump power in watts; 
wp is the pump beam size; 

1 is the loss coefficient on the pump wavelength (in 1/cm); 
loss is the out-put loss coefficient; 
am is the maximum gain coefficient; 

Per is the critical power of the self-focusing in the active medium 
(in watts); 

Dam is the group- velocity dispersion of the active medium (in fs~2/cm, 
normal dispersion has a negative sign); 
Is is the gain saturation intensity; 

q is the complex Gaussian beam parameter; 
Pw is the pulse power; 

delta is the pulse width (for Gaussian pulse); 
psi is the pulse chirp; 

alpha and beta is the generation beam parameters (see description) 
Lengths are given in centimeters; powers are given in watts 

COMPLEX*16 i,qs,q(15),qq,t ! i is the imaginary unit, qs is the initial 

REAL*8 a,c,bl,b2,b,f,z,epsl,eps2,ro(15),n,x,ka(4),qu(4),tau(4) 
REAL*8 delta(4),psi(4),del,delta0,ps,psi0,tg 

REAL*8 lambda, Pi, bis, bs,bf,betaO,alphaO,k, beta, alpha,kappa,wp,P,l 
REAL*8 gamma, gammaO, loss, Pw,g, Is, S, am, Per, Pd(l),Der, Dam, Lcav,Tcav 
INTEGER Num 

DATA a,c,z,f,n,epsl/30d0,60d0,0.28d0,15d0,2.442d0,ld-3/ 

DATA bls,bs,bf,wp,P,l,loss/0d0,20d0,55d0,100d-4,1.5d0,6d0,5d-2/ 

DATA Pi,lambda/3.14159265358979323846264338d0,2.5d-4/ 

DATA am,Pcr,cv,Dam,Is/9dO,965d3,3dlO,-2054dO,ld4/ 

i=dcmplx(0.,l.) 

open(l,file='ro.dat') 

open(2,filc='der.dat') 

write(l,*)a,c,wp,P 

! w0=100 mkm is the initial size of the plane wave: 
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qs = (0.,. 7957747152) 

! wmax=10 cm is maximum size of the simulated mode 

eps2 = ld2*Pi/lambda 

k = 2dO*Pi/lambda ! is the wave number 

! It takes into account the pump wave damping in the active medium: 
kappa =.5506035739d0*P 

S = 2d0*0.48d0*P*.2900302115d-4/Pi/wp**2 ! is the pump parameter 
! sigma.a*T_r/h/nu_p=.2900302115d-4. 
! 0.48 - averaging along propagation axis. 
! sigma_a, T_r, nu_p are the absorption 
! cross section of the active medium, 
! the gain relaxation time, and the gain 
! wavelength, respectively 

x = z/n ! is the optical length of the gain medium 

tg = 4d0/z ! is the gain band width in fs/cm 

dx = x/ld3 ! is the step size 

DO 11=1,201 ! scanning on b 

b = bs + (bf-bs)*(Il-l)/200dO 
blf = b-z 

DO 12=1,201 ! scanning on bl 

bl = bis + (blf-bls)*(I2-l)/200dO 
b2 = b-bl-z 

Lcav = a + b + c ! is the cavity length 

Tcav = (2dO*Lcav/cv)*ldl5 ! is the cavity period fs 

write(*,*)Il,I2 ! just step numbers ! 

gammaO = OdO 

Pw = lOdO ! initial pulse power in watts 

deltaO = ld3 ! initial pulse width in fs 

psiO = OdO ! initial pulse chirp 

c beam propagation 

Num = 

qq = qs 
2 continue 

Num = Num + 1 

Pwold = Pw 

c ABCD-modul: output mirror - folding mirror - active medium 
q(l) = qq 
q(2) = q(l) + a 
q(3) = Id0/(ld0/q(2) - IdO/f) 
q(4) = q(3) + bl 



if(-ld0/dimag(ld0/q(4)).lt.0d0)goto 1 

betaO = -k*dreal(ld0/q(4))/2d0 

alphaO = dsqrt(-ld0/dimag(ld0/q(4))/k) 
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g = am*S/(ldO + S + 
#Pw*dcltaO*dsqrt(Pi/2dO)/(Pi*alphaO**2*Is*Tcav)) ! saturated gain 

c Active medium (Runge-Kutta fourth-order method) 

DO J=l,1000 

ka(l) = (-0.5d0/alpha0**4 + 2d0*beta0**2 + 
#2dO*k*kappa*dexp(-l*dx*(J-l))/wp**2)/k + 

#(dsqrt(2dO)/Pi)*(Pw/Pcr)*dexp(2*gamniaO)/alphaO**4/k !beta 

qu(l) = - 2d0*beta0*alpha0/k - 
#2d0*g*alpha0**3/(ld0+2d0*alpha0**2/wp**2)**(3/2)/wp' ,! *2 lalpha 

tau(l) = g/(ldO + 2d0*alpha0**2/wp**2) - 2d0*tg**2/delta0**2 + 

! gamma 

#Dam*psiO 

delta(l) = 2d0*tg**2*(ld0/delta0- delta0**3*psi0**2) - ! delta 

#2dO*Dam*deltaO*psiO 

psi(l) = 2dO*Dam*(psiO**2 - Id0/delta0**4) - 8d0*tg**2*psi0/ !psi 
#delta0**2 + (2dO/Pi)*(Pw/Pcr)*dexp(2*gammaO)/alphaO**2/deltaO**2 



ka(2) = (-0.5d0/(alpha0+0.5d0*dx*qu(l))**4 + 
#2d0*(beta0+0.5d0*dx*ka(l))**2 + 
#2dO*k*kappa*dexp(d*dx*(J-0.5dO))/wp**2)/k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)*dexp(2*gammaO+0.5dO*dx*tau(l))/ 
#(alpha0+0.5d0*dx*qu(l))**4/k 

qu(2) = - 2d0*(beta0+0.5d0*dx*ka(l))*(alpha0+0.5d0*dx*qu(l))/k- 
#2d0*g*(alpha0+0.5d0*dx*qu(l))**3/(ld0+ 
#2d0*(alpha0+0.5d0*dx*qu(l))* s,! 2/wp s,! *2)**(3/2)/wp s,!S,! 2 

tau(2) = g/(ldO + 2d0*(alpha0+0.5d0*dx*qu(l))**2/wp**2) - 
#2dO*tg**2/(deltaO+0.5dO*dx*delta(l))**2+Dam*(psiO+0.5dO*dx*psi(l)) 

delta(2) = 2d0*tg**2*(ld0/(dclta0+0.5d0*dx*delta(l)) - 
#(deltaO+0.5dOMxMelta(l))**3*(psiO+0.5dO*dx*psi(l))**2) - 
#2dO*Dam* (dclta0+0.5d0*dx*dclta(l ) ) * (psi0+0.5d0*dx*psi(l ) ) 

psi(2) = 2dO*Dam*((psiO+0.5dO*dx*psi(l))**2- 
#ldO/ (delta0+0.5d0*dx*delta(l ) )**4) - 
#8d0*tg**2*(psi0+0.5d0*dx*psi(l))/ 
#(delta0+0.5d0*dx*delta(l))**2 + (2dO/Pi)*(Pw/Pcr)* 
#dcxp(2*(gammaO+0.5dO*dx s,! tau(l)))/(alphaO+0.5dO s,! dx s,! qu(l))**2/ 
#(delta0+0.5d0*dx*delta(l))**2 

ka(3) = (-0.5d0/(alpha0+0.5d0*dx*qu(2))**4 + 
#2d0*(beta0+0.5d0*dx*ka(2))**2 + 
#2dO*k*kappa*dcxp(d*dx*(J-0.5dO))/wp**2)/k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)Mexp(2*gammaO+0.5dO*dx*tau(2))/ 
#(alpha0+0.5d0*dx*qu(2))**4/k 

qu(3) = - 2dO*(betaO+0.5dOMx*ka(2))*(alphaO+0.5dO*dx*qu(2))/k- 
#2d0*g*(alpha0+0.5d0*dx*qu(2))**3/(ld0+ 
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#2d0*(alpha0+0.5d0*dx*qu(2))**2/wp**2)**(3/2)/wp**2 

tau(3) = g/(ldO + 2d0*(alpha0+0.5d0*dx*qu(2))**2/wp**2) - 

#2d0*tg**2/(delta0+0.5d0*dx*delta(2))**2+Dam*(psi0+0.5d0*dx*psi(2)) 
delta(3) = 2d0*tg**2*(ld0/(delta0+0.5d0*dx*delta(2)) - 

#(dclta0+0.5d0*dx*delta(2))**3*(psi0+0.5d0*dx*psi(2))**2) - 

#2dO*Dam*(deltaO+0.5dO*dx*delta(2))*(psiO+0.5dO*dx*psi(2)) 

psi(3) = 2dO*Dam*((psiO+0.5dO*dx*psi( 2 ))** 2 - 
#ld0/(delta0+0.5d0*dx*delta(2))**4) - 
#8d0*tg**2*(psi0+0.5d0*dx*psi(2))/ 
#(delta0+0.5d0*dx*delta(2))**2 + (2dO/Pi)*(Pw/Pcr)* 
#dcxp(2*(gamma0+0.5d0*dx*tau(2)))/(alpha0+0.5d0 ,,! dx*qu(2))**2/ 
#(delta0+0.5d0*dx*delta(2))**2 

ka(4) = (-0.5d0/(alpha0+dx*qu(3))**4 + 
#2d0*(beta0+dx*ka(3))**2 + 
#2d0*k*kappa*dexp(-l*dx* J) /wp**2) /k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)*dexp(2*gammaO+dx*tau(3))/ 
#(alpha0+dx*qu(3))**4/k 

qu(4) = - 2d0*(beta0+dx*ka(3))*(alpha0+dx*qu(3))/k- 
#2d0*g*(alpha0+dx*qu(3))**3/(ld0+ 
#2d0*(alpha0+dx*qu(3))**2/wp**2)**(3/2)/wp**2 

tau(4) = g/(ldO + 2d0*(alpha0+dx*qu(3))**2/wp**2) - 
#2d0*tg**2/(delta0+dx*delta(3))**2 + Dam*(psi0+dx*psi(3)) 

delta(4) = 2d0*tg**2*(ld0/(delta0+dx*delta(3)) - 
#(delta0+dx*delta(3))**3*(psi0+dx*psi(3))**2) - 
#2dO*Dam*(dcltaO+dx*delta(3))*(psiO+dx*psi(3)) 

psi(4) = 2dO*Dam*((psiO+dx*psi(3))**2- 
#ld0/(delta0+dx*delta(3))**4) - 
#8d0*tg**2*(psi0+dx*psi(3))/ 
#(delta0+dx*delta(3))**2 + (2dO/Pi)*(Pw/Pcr)* 
#dcxp(2*(gamma0+dx*tau(3)))/(alpha0+dx*qu(3))**2/ 
#(delta0+dx*delta(3))**2 

beta = betaO + (dx/6d0)*(ka(l)+2d0*ka(2)+2d0*ka(3)+ka(4)) 

alpha = alphaO + (dx/6d0)*(qu(l)+2d0*qu(2)+2d0*qu(3)+qu(4)) 

gamma = gammaO + (dx/6d0)*(tau(l)+2d0*tau(2)+2d0*tau(3)+tau(4)) 

dcl=dclta0+(dx/6d0)*(delta(l)+2d0*delta(2)+2d0*delta(3)+delta(4)) 

P s=psi0+(dx/6d0)*(psi(l)+2d0*psi(2)+2d0*psi(3)+psi(4)) 

if (alpha. It. 0. or. alpha. gt. 10. or. gamma.lt. 0. or. gamma. gt. 10. 
#or.del.lt.0.)goto 1 

alphaO = alpha 

betaO = beta 

gammaO = gamma 

deltaO = del 

psiO = ps 

END DO 

Pw = Pw*dexp(2d0*gamma0) 
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gammaO = OdO 
if(Pw.le.ld-10.or.Pw.gt.ld20)goto 1 
t = k/(-2d0*beta0 - i/alpha0**2) 
q(5) = dcmplx(dreal(t),dimag(t)) 
c 

c ABCD-modul: active medium - second folding mirror - second flat mirror 
(and backwards) 

q(6) = q(5) + b2 

q(7) = Id0/(ld0/q(6) - IdO/f) 

q(8) = q(7) + c 

q(9) = q(8) + c 

q(10) = Id0/(ld0/q(9) - IdO/f) 

q(ll) = q(10) + b2 



if(-ldO/dimag(ldO/q(ll)).lt.OdO)goto 1 

betaO = -k*dreal(ld0/q(ll))/2d0 

alphaO = dsqrt(-ldO/dimag(ldO/q(ll))/k) 

g = am*S/(ldO + S + 
#Pw*dcltaO*dsqrt(Pi/2dO)/(Pi*alphaO**2*Is ,,! Tcav)) 

c Active medium (Runge-Kutta fourth-order method) 
DO J=l,1000 

ka(l) = (-0.5d0/alpha0**4 + 2d0*beta0**2 + 
#2dO*k*kappa*dexp(-l*dx*(J-l))/wp**2)/k + 
# (dsqrt (2d0) /Pi) * (Pw/Pcr) *dexp(2*gamma0) /alphaO* *4/k 

qu(l) = - 2d0*beta0*alpha0/k - 
#2d0*g*alpha0**3/(ld0+2d0*alpha0 s,! *2/wp s,! *2) s,! *(3/2)/wp s,!S,! 2 

tau(l) = g/(ldO + 2d0*alpha0**2/wp**2) - 2d0*tg**2/delta0**2 + 
#Dam*psiO 

delta(l) = 2d0*tg**2*(ld0/delta0- delta0**3 ,,! psi0**2) - 
#2dO*Dam*deltaO*psiO 

psi(l) = 2dO*Dam*(psiO**2 - Id0/delta0**4) - 8d0*tg**2*psi0/ 
#delta0**2 + (2dO/Pi)*(Pw/Pcr)*dexp(2*gammaO)/alphaO**2/deltaO**2 



ka(2) = (-0.5d0/(alpha0+0.5d0*dx*qu(l))**4 + 
#2d0*(beta0+0.5d0*dx*ka(l))**2 + 
#2dO*k*kappa*dexp(-l*dx*(J-0.5dO))/wp**2)/k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)*dexp(2*gammaO+0.5dO*dx*tau(l))/ 
#(alpha0+0.5d0*dx*qu(l)) s,! *4/k 

qu(2) = - 2d0*(beta0+0.5d0*dx*ka(l))*(alpha0+0.5d0*dx*qu(l))/k- 
#2d0*g*(alpha0+0.5d0*dx*qu(l))**3/(ld0+ 
#2d0*(alpha0+0.5d0*dx*qu(l))**2/wp**2)**(3/2)/wp**2 

tau(2) = g/(ldO + 2d0*(alpha0+0.5d0*dx*qu(l))**2/wp**2) - 
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#2dO*tg**2/(deltaO+0.5dO*dx*delta(l))**2+Dam*(psiO+0.5dO*dx*psi(l)) 

delta(2) = 2d0*tg**2*(ld0/(delta0+0.5d0*dx*delta(l)) - 
#(delta0+0.5d0*dx*delta(l))**3*(psi0+0.5d0*dx*psi(l))**2) - 
#2dO*Dam* (dclta0+0.5d0*dx*delta(l ) ) * (psi0+0.5d0*dx*psi(l ) ) 

psi(2) = 2dO*Dam*((psiO+0.5dO*dx*psi(l))**2- 
#ldO/ (deltaO+0. 5d0*dx*delta( 1 ) ) **4) - 
#8d0*tg**2*(psi0+0.5d0*dx*psi(l))/ 
#(delta0+0.5d0*dx*delta(l))**2 + (2dO/Pi)*(Pw/Pcr)* 
#dcxp(2*(gamma0+0.5d0*dx*tau(l)))/(alpha0+0.5d0 ,,! dx*qu(l))**2/ 
#(delta0+0.5d0*dx*delta(l))**2 

ka(3) = (-0.5d0/(alpha0+0.5d0*dx*qu(2))**4 + 
#2d0*(beta0+0.5d0*dx*ka(2))**2 + 
#2dO*k*kappa*dcxp(-l*dx*(J-0.5dO))/wp**2)/k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)Mexp(2*gammaO+0.5dO*dx*tau(2))/ 
#(alpha0+0.5d0*dx*qu(2))**4/k 

qu(3) = - 2d0*(beta0+0.5d0*dx*ka(2))*(alpha0+0.5d0*dx*qu(2))/k- 
#2d0*g*(alpha0+0.5d0*dx*qu(2))**3/(ld0+ 
#2d0*(alpha0+0.5d0*dx*qu(2))**2/wp**2)**(3/2)/wp**2 

tau(3) = g/(ldO + 2d0*(alpha0+0.5d0*dx*qu(2))**2/wp**2) - 
#2d0*tg**2/(delta0+0.5d0*dx*delta(2))**2+Dam*(psi0+0.5d0*dx*psi(2)) 

delta(3) = 2d0*tg**2*(ld0/(dclta0+0.5d0*dx*delta(2)) - 
#(dclta0+0.5d0*dx*delta(2))**3*(psi0+0.5d0*dx*psi(2))**2) - 
#2dO*Dam*(deltaO+0.5dOMx*delta(2))*(psiO+0.5dO*dx*psi(2)) 

psi(3) = 2dO*Dam*((psiO+0.5dO*dx*psi(2))**2- 
#ld0/(delta0+0.5d0*dx*delta(2))**4) - 
#8d0*tg**2*(psi0+0.5d0*dx*psi(2))/ 
#(delta0+0.5d0*dx*delta(2))**2 + (2dO/Pi)*(Pw/Pcr)* 
#dcxp(2*(gamma0+0.5d0*dx*tau(2)))/(alpha0+0.5d0 s,! dx s,! qu(2))**2/ 
#(delta0+0.5d0*dx*delta(2))**2 

ka(4) = (-0.5d0/(alpha0+dx*qu(3))**4 + 
#2d0*(beta0+dx*ka(3))**2 + 
#2d0*k*kappa*dexp(-l*dx* J) /wp**2) /k + 
#(dsqrt(2dO)/Pi)*(Pw/Pcr)*dexp(2*gammaO+dx*tau(3))/ 
# (alpha0+dx*qu(3) ) * *4/k 

qu(4) = - 2d0*(beta0+dx*ka(3))*(alpha0+dx*qu(3))/k- 
#2d0*g*(alpha0+dx*qu(3))**3/(ld0+ 
#2d0*(alpha0+dx*qu(3))**2/wp**2)**(3/2)/wp**2 

tau(4) = g/(ldO + 2d0*(alpha0+dx*qu(3))**2/wp**2) - 
#2d0*tg**2/(delta0+dx*delta(3))**2 + Dam*(psi0+dx*psi(3)) 

delta(4) = 2d0*tg**2*(ld0/(delta0+dx*delta(3)) - 
#(delta0+dx*delta(3))**3*(psi0+dx*psi(3))**2) - 
#2dO*Dam*(deltaO+dx*delta(3))*(psiO+dx*psi(3)) 

psi(4) = 2dO*Dam*((psiO+dx*psi(3))**2- 
#ld0/(delta0+dx*delta(3))**4) - 
#8d0*tg**2*(psi0+dx*psi(3))/ 
#(delta0+dx*delta(3))**2 + (2dO/Pi)*(Pw/Pcr)* 
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#dexp(2*(gamma0+dx*tau(3)))/(alpha0+dx*qu(3))**2/ 

#(delta0+dx*delta(3))**2 

beta = betaO + (dx/6d0)*(ka(l)+2d0*ka(2)+2d0*ka(3)+ka(4)) 
alpha = alphaO + (dx/6d0)*(qu(l)+2d0*qu(2)+2d0*qu(3)+qu(4)) 
gamma = gammaO + (dx/6d0)*(tau(l)+2d0*tau(2)+2d0*tau(3)+tau(4)) 
del=delta0+(dx/6d0)*(delta(l)+2d0*delta(2)+2d0*delta(3)+delta(4)) 
P s=psi0+(dx/6d0)*(psi(l)+2d0*psi(2)+2d0*psi(3)+psi(4)) 
if (alpha. It. 0. or. alpha, gt. 10. or. gamma.lt. 0. or. gamma.gt. 10. 

#or.del.lt.0.)goto 1 
alphaO = alpha 
betaO = beta 
gammaO = gamma 
deltaO = del 
psiO = ps 

END DO 



Pw = Pw*dcxp(2d0*gamma0) 

gammaO = OdO 
if(Pw.lc.ld-10.or.Pw.gt.ld20)goto 1 
t = k/(-2d0*beta0 - i/alpha0**2) 
q(12) = dcmplx(drcal(t),dimag(t)) 
c ABCD-modul for the residuary propagation up to out-put mirror 
q(13) = q(12) + bl 
q(14) = Id0/(ld0/q(13) - IdO/f) 
q(15) = q(14) + a 
Pw = Pw*dexp(-loss) 



if(Pw.le.ld-10.or.Pw.gt.ld20)goto 1 ! criteria for the power 
if(Num.gt.5000)goto 1 ! and iteration number 

c w"2*Pi/lambda converts the initial part of the beam parameter to the 
beam size 

DO 14=1,15 

ro(I4) = -Id0/dimag(ld0/q(14)) 
END DO 

c Beam is to have the positive size and hasn't to be too large 
if(ro(l).le.0..or.ro(2).le.0..or.ro(3).le.0..or. 

2 ro(4).le.0..or.ro(5).le.0..or.ro(6).le.0..or. 

3 ro(7).le.0..or.ro(8).lc.0..or.ro(9).le.0..or. 

4 ro(10).lc.0..or.ro(ll).le.0..or.ro(12).le.0..or. 

5 ro(13).le.0..or.ro(14).lc.0..or.ro(15).le.0.)goto 1 
if(ro(l).gt.eps2.or.ro(2).gt.eps2.or.ro(3).gt.eps2.or. 

2 ro(4).gt.cps2.or.ro(5).gt.eps2.or.ro(6).gt.eps2.or. 

3 ro(7).gt.eps2.or.ro(8).gt.eps2.or.ro(9).gt.eps2.or. 
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4 ro(10).gt.eps2.or.ro(ll).gt.eps2.or.ro(12).gt.eps2.or. 

5 ro(13).gt.eps2.or.ro(14).gt.eps2.or.ro(15).gt.eps2)goto 1 
qq = q(15) 

c Pulse power stability 

if(abs(ro(15)/ro(l)-ld0).gt.epsl. 
#or.abs(Pw/Pwold-ldO).gt.epsl)goto 2 
c Out-put for the stable pulse 

write(l,*)b,bl,Pw,sqrt(ro(5)*lambda/Pi),deltaO,Num 
1 continue 

END DO 
END DO 

close(l) 
close(2) 
STOP 
END 
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